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1.  Introduction 

Our  research  project  with  ARO  started  with  a  small  effort  in  1977  to  explore 
certain  connections  we  have  noted  between  inverse-scattering  theory  (Gelfand-Levitan, 
Marchenko,  etc.)  and  linear  least-squares  estimation.  As  indicated  by  our  list  of 
publications  with  ARO  support,  this  has  been  a  very  fruitful  area  of  research.  We  shall 
not  summarize  all  our  results  here;  we  note,  however,  that  one  important  result  is  our 
establishment  of  a  fundamental  correspondence  between  one-dimensional  inverse 
scattering  in  a  lossless  medium  and  efficient  triangular  factorization  of  certain 
structured  (close-to-Toeplitz)  matrices.  Then  by  broadening  the  notion  of  structured 
matrices  we  were  able  to  recast  many  signal  processing  problems  as  (non-classical) 
inverse  scattering  problems,  e.g.,  partial  realization,  decoding  of  error-correcting  codes, 
and  polynomial  zero-location.  This  led  us  to  develop  a  general  scheme  (see  Lev-Ari 
and  Kailath  (1986))  that  can  be  applied  to  solve  each  one  of  these  generalized  inverse 
scattering  problems  and,  at  the  same  time,  to  compute  the  triangular  factors  of  a 
certain  associated  structured  matrix. 

Our  work  has  led  us  to  realize  that  these  seemingly  unrelated  problems  can  be 
grouped  into  two  distinct  families: 

(1)  Inverse  scattering  problems,  which  involve  the  reconstruction  of  internal 
parameters  of  an  unknown  system  from  given  boundary  (i.e.,  input/output)  data. 

(2)  Matrix  factorization  problems,  which  involve  obtaining  the  inertia  and  often 
also  the  triangular  (lower-diagonal-upper)  factorization  of  a  given  Hermitian 
matrix. 

We  showed  that  there  exists  a  direct  correspondence  oe'ween  these  two  families 
of  problems  via  the  notion  of  energy  conservation  or  losslessness.  For  this,  we  first 
introduced  a  procedure  for  associating  Hermitian  matrices  with  incident  and  reflected 
waves  in  a  lossless  cascade  network  and  then  showed  that  the  information  required  to 
construct  the  triangular  factorization  of  these  matrices  was  provided  by  the  internal 
signals  in  the  network.  On  the  other  hand,  we  focused  on  a  particular  procedure  for 
solving  the  inverse  scattering  problem,  known  as  ‘layer  peeling’  [Bruckstein  and 
Kailath  (1987)].  Here,  the  unknown  system  is  identified  by  layers  or  sections:  the 
outermost  layer  is  identified  first,  then  peeled  off  to  reveal  the  input/output  signals 
associated  with  the  next  layer.  Translating  this  procedure  to  matrix  language  results  in 
an  extension  of  the  well-known  Cholesky  factorization  method:  the  triangular  factor  of 
a  given  matrix  is  computed  recursively,  column-by-column.  Each  time  a  new  column 
has  been  computed  its  effect  is  removed  by  subtracting  a  suitable  matrix  of  rank  one. 
The  recursive  nature  of  these  layer-peeling  and  factorization  procedures  results  in  a 
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cascade  structure  for  the  associated  network  models.  Thus  both  families  of  problems 
are,  in  fact,  two  facets  of  the  same  phenomenon,  and  the  same  relation  exists  between 
the  recursive  procedures  used  to  solve  these  problems. 

Moreover,  with  an  appropriate  formulation  we  were  able  to  extend  our  conceptual 
framework  to  cascade  networks  that  display  a  rather  generalized  concept  of 
losslessness  [Lev-Ari  and  Kailath  (1986)].  It  may  be  useful  to  recall  here  that  cascade 
network  models  abound  in  the  circuits  and  systems  literature.  They  are  central  to  the 
classical  filter  synthesis  methods  of  Darlington  and  Brune,  as  well  as  to  the  modem 
filter  synthesis  methods  of  Deprettere  and  Dewilde  (1980),  Rao  and  Kailath  (1984), 
and  Vaidyanathan  and  Mitra  (1984).  They  arise  in  the  stochastic  filtering  of  stationary 
and  nonstationary  processes,  as  well  as  in  the  related  problems  of  efficient  factorization 
and  inversion  of  Toeplitz  and  near-Toeplitz  matrices  [see,  e.g.,  Lev-Ari  and  Kailath 
(1984),  and  Kailath  (1986)]  and  extensions  thereof  [Lev-Ari  and  Kailath  (1986)].  They 
are  implicit  in  much  of  the  research  on  fast  algorithms  for  Toeplitz  and  Hankel 
matrices  [see,  e.g.,  Friedlander,  Morf,  Kailath  and  Ljung  (1979),  and  Heinig  and  Rost 
(1984)].  Since  the  layered  earth  model  of  Gopilaud  is  a  lossless  cascade  network,  the 
inverse  scattering  literature  is  also  concerned  with  such  networks  [see  e.g.,  Robinson 
(1975),  Bruckstein  and  Kailath  (1987),  Brucksteirv,  L6vy  and  Kailath  (1985)].  Finally, 
such  networks  arise  in  the  solution  of  several  problems  in  system  theory,  namely  in  the 
construction  of  efficient  tests  for  stability  and  zero  location  [see,  e.g..  Jury  (1964)],  in 
the  solution  of  the  partial  realization  problem  [Citron,  Bruckstein  and  Kailath  (1984)] 
and  of  the  stochastic  realization  problem  [Mullis  and  Roberts  (1986)],  and  recently 
also  in  model  order  reduction  techniques  [Ball  and  Gohberg  (1986)]. 


Ifl 

w2 
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Vi 

V2 

Fig.  1  Lossless  Cascade  Network 


It  will  also  be  helpful  to  be  more  specific  about  the  connections  between  cascade 
networks,  inverse-scattering  and  matrix  factorization.  Consider  a  cascade  of  sections 
(or  layers),  with  p  discrete-time  inbound  signals  u  ,  flowing  from  the  input 
terminals  into  the  system,  and  q  output  signals  v  flowing  out  of  the  system  to  the 
output  terminals  (Fig.  1).  This  cascade  system  is  called  (internally)  lossless  if  the 
energy  (=  /  2  norm)  input  into  each  section  equals  the  energy  output  from  the  same 
section,  viz., 


I|u«-ill2  +  IM2  =  I|ut  H2  -f-  Hv^ill2 


i  >  1 


where  u,  ,  v(-  denote  the  signals  flowing  through  the  boundary  between  the  i  -th  and 
the  (z'-f-l)-th  layer.  Alternatively,  we  may  characterize  losslessness  in  terms  of  the 
scattering  matrices  Z;(z),  which  relate  the  Z-transforms  of  these  signals,  i.e., 

*  r  -v 

M,(z)  M;_  i(z) 

U:1(2)J  -*«[  ;i(z)j  <>-2> 

It  is  well-known  that  a  necessary  and  sufficient  condition  for  losslessness  is 

It(z)|X(l/z*)]*  =  /  (1.3) 

where  the  asterisk  (*)  denotes  Hermitian  transpose  (complex  conjugate  for  scalars). 
Often  we  need  to  use  the  chain-scattering  matrices  0,  (z ),  which  relate  the  signals  at 
one  side  of  each  layer  to  the  signals  on  the  other  side,  i.e.. 


v,.(z)  -  W‘(z)  v.^z) 


i  £  1 


Since  the  losslessness  constraint  can  be  rewritten  in  the  form 
llu.  ll2  -  ||v,||2  =  ||u,-_j||2  —  ||v,-_1||2  it  follows  that  the  chain-scattering  matrices  are 
J-orthogonal,  i.e., 


0,(z)/[0(.(l/z*)]  =J  ,  J  :=  diag{lp  ,-Iq}  (1.5) 

The  layer  peeling  procedure  can  now  be  described  as  follows: 

Beginning  with  the  given  boundary  data  u0(z)  and  v0(z)  repeat  the 
following  steps  for  /  =  1,  2,  .  .  . 

(a)  Identify  the  internal  parameters  of  the  /-th  section  from  u^^z) , 

v.-iOO. 

(b)  Compute  m,  (z),  v,  (z)  via  (1.4). 

Notice  that  the  losslessness  property  is  not  required  in  order  to  carry  out  this 
procedure,  only  that  we  are  able  to  identify  each  section  from  its  input/output  data. 
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This  observation,  first  made  by  Bruckstein  and  Kailath  (1987),  indicates  that  the  layer¬ 
peeling  procedure  may  apply  to  cascade  networks  that  do  not  satisfy  the  standard 
notion  of  losslessness  (1.3). 


Losslessness  is,  nevertheless,  essential  to  the  association  of  inverse  scattering 
problems  with  matrix  factorization  problems,  as  demonstrated  by  Kailath,  Bruckstein 
and  Morgan  (1986).  They  analyze  the  layer- peeling  procedure  associated  with  a 
lossless  transmission-line  model,  where 


®i(z) 


1 

V  i  - 1  ki  i2 


(1.6) 


is  a  chain-scattering  matrix  with  a  single  storage  element  (z).  They  show  that  the 
principle  of  energy  conservation  results  in  an  identity  of  the  form 

L(«o)L  (mq)  —  L(v0)L  (vq)  =  XX  (1-7) 


where  X  is  a  lower-triangular  matrix  whose  y'-th  column  consists  of  the  samples  of 
the  signal  stored  in  the  y'-th  section  of  the  network,  and  L(m0)  (resp.  L(v0))  is  a 
lower-triangular  Toeplitz  matrix  constructed  from  the  coefficients  of  the  Z -transform 
u0(z)  (resp.  v0(z)),  namely 


L(0 


f  0 

f  \  f  0  O 

fl  f  1  /  0 


(1.8a) 


where 


/(O  :=  f/i*1' 


i=0 


(1.8b) 


is  a  Z -transform  of  a  signal  f  ={/,•;  0  <  i  <  <»}.  Thus,  the  Hermitian  matrix 


R  :=  L(h0)L*(u0)  -  L(v0)L*(v0)  (1.9) 


can  be  factored  by  applying  an  appropriate  layer-peeling  procedure  to  the  pair  of 
signals  {«o(z)»  vo(z)}  and  reading  off  the  elements  of  the  triangular  factor  X  from 
the  signals  stored  in  the  delay  element  in  each  section. 

The  association  between  inverse-scattering  problems  and  efficient  matrix 
factorization  has  motivated  us  to  look  for  a  similar  association  involving  structured 
matrices  that  are  not  of  the  form  (1.9),  e.g.,  Hankel  matrices,  sums  of  Toeplitz  and 
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Hankel  matrices,  Bezoutians  with  respect  to  various  curves,  and  Vandermonde 
matrices.  The  key  to  the  extension  of  the  results  (1.7)-(1.9)  beyond  the  standard  notion 
of  losslessness  is  the  observation  that  the  matrix  R  in  (1.9)  has  a  very  specific 
displacement  structure  [see,  e.g.,  Lev-Ari  and  Kailath  (1984),  and  Kailath  (1986)]. 
One  convenient  way  to  describe  this  structure  is  via  the  bivariate  Z-transform  (= 
generating  function)  of  R,  viz., 

R(z,w )  :=  [l  z  z2  .  .  .  ]  R  [  1  w  w2...]*  (1.10) 


Applying  this  definition  to  (1.9)  we  obtain 


R(z,w) 


“o(*W  (w)  -  Vq(z)Vq  (w) 
1  -  zw* 


(1.11) 


which,  for  various  choices  of  {u0(z),  v0(z)}  describes  the  so-called  quasi-Toeplitz 
family  of  matrices  [see  Lev-Ari  and  Kailath  (1986)].  The  quasi-Hankel  family  of 
matrices,  which  arises  in  the  study  of  Bezoutians,  also  has  a  simple  form,  viz., 

R  (z  ,w )  =  .gfrAVl-jjUgV)  (,.12) 

z  —  w 


A  comparison  of  (1.11)  with  (1.12),  as  well  as  our  previous  work  in  this  area  suggests 
a  generalization  to  matrices  whose  generating  functions  are  of  the  form 


R(z,w) 


G(z)J  GV) 
d{z  ,w ) 


(1.13) 


where  G(z )  is  a  row  vector  of  functions,  J  is  a  constant  nonsingular  Hermitian 
matrix  of  suitable  dimensions  and  the  displacement  function  d(z,w)  is  a  fixed 
function  with  the  Hermitian  property  d*(z  ,w)  =  d  (w  ,z ).  The  triplet 
{d(z  ,w ),  G  (z ),  J  }  has  been  called  a  generator  of  the  Hermitian  matrix  R,  since  it 
uniquely  determines  R(z  ,w)  (and  hence  also  R)  via  (1.13)  [for  more  details,  see 
Lev-Ari  and  Kailath  (1986)]. 


In  Section  2  we  review  the  generalized  factorization  algorithm  we  developed  for 
structured  matrices  with  generating  functions  of  the  form  (1.13).  However,  these 
results  were  developed  for  strongly-regular  matrices,  viz.,  those  with  non-vanishing 
leading  singular  minors.  Such  singularities  have  not  been  easy  to  handle,  even  in  the 
special  cases  of  Toeplitz  and  Hankel  matrices,  though  much  more  progress  has  been 
made  for  Hankel  matrices.  Our  generating  function  approach  seems  to  allow  the 
development  of  efficient  procedures,  which  involve  only  scalar  computations,  for 
overcoming  singularities  in  the  factorization.  In  contrast,  the  occurrence  of 
singularities  in  non-structured  matrices  can  be  overcome  only  by  inverting  a  matrix 
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whose  size  is  proportional  to  the  depth  of  the  singularity. 

In  another  direction,  we  have  noticed  that  our  results  can  be  applied  to  the  much- 
studied  (for  nearly  150  years)  problems  of  determining  the  stability  or,  more  generally, 
the  root-distribution  with  respect  to  a  given  curve  in  the  complex  plane  (e.g.,  the 
imaginary  axis  or  the  unit  circle).  This  is  because  by  a  classical  result  of  Hermite 
(1856),  the  root-  distribution  is  determined  by  the  inertia  of  a  certain  structured  matrix 
(known  as  Bezoutian),  which  has  displacement  structure  in  our  sense  (see  (1.11)- 

(1.13) ).  It  turns  out  that  our  general  factorization  scheme  readily  yields  the  classical 
Routh-Hurwitz  procedure  for  the  imaginary  axis  and  the  classical  Schur-Cohn 
procedure  for  the  unit  circle  [see  Lev-Ari,  Bistritz  and  Kailath  (1987)].  Moreover,  as 
will  be  described  in  more  detail  in  Section  3,  our  approach  shows  that  there  exist 
many  alternatives  to  these  classical  procedures,  some  of  which  are  more  efficient  than 
the  classical  ones.  Our  work  provides  an  organized  and  complete  approach  to  this 
much-studied  subject,  relating  old  tests  to  each  other  and  indicating  where  new  tests 
remain  to  be  developed  and  explored. 

Finally,  in  Section  4  we  describe  yet  another  direction  of  research  on 
displacement  structure  that  leads  to  several  new  results,  especially  on  orthogonal  (QR) 
rather  than  triangular  (LDU)  factorization  and  correspondingly  on  least-squares  rather 
than  exact  solutions  of  linear  equations.  The  main  point  here  is  to  go  back  to  the 
matrix  form  of  the  generating  function  representation  (1.13),  by  ‘inverting’  the  Z- 
transform  operation  (1.10).  For  instance,  the  resulting  matrix  equivalent  of  (1.11)  is 

R  -  Z  R  Z*  =  UqUq  -  v0Vq  (1.14) 

where  Z  is  the  lower  shift  matrix  with  l’s  on  the  first  subdiagonal  and  0’s 
elsewhere,  and  u0  ,  v0  are  column  vectors  consisting  of  the  coefficients  of  the  power 
series  expansions  of  u0  ,  v0  ,  respectively.  It  turns  out  that  the  main  features  of  our 
structured  factorization  procedure  are  preserved  even  when  we  replace  Z  and  Z*  in 

(1.14)  by  arbitrary  (lower-triangular)  matrices  Fj  and  F2  ,  so  that 

R-F^F^  =  uu*  -  vv*  (1.15) 

A  matrix  R  that  satisfies  this  equation  admits  a  structured  factorization  procedure 
even  though  (1.15)  does  not  correspond  any  more  to  a  simple  generating  function 
representation  like  (1.13).  In  fact,  a  general  (efficient)  factorization  scheme  can  be 
constructed  for  all  matrices  R  that  satisfy  a  displacement  equation  of  the  form 

Fj  R  F2  —  F3  RF4  =  uu*  —  vv*  (1.16) 

where  F,  for  1  =1,2, 3, 4  are  arbitrary  lower-triangular  matrices.  Furthermore,  the 
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same  general  scheme  can  be  used  to  obtain  efficient  algorithms  for  matrix  inversion 
(exact  and  least-squares)  and  QR  factorization;  this  is  achieved  by  applying  the  general 
factorization  procedure  to  a  suitable  composite  matrix,  which  is  constructed  from  the 
given  matrix  R.  We  have  recently  begun  to  explore  the  application  of  such  composite 
matrix  formulations  to  new  problems,  as  described  in  Section  4. 
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2:  Factorization  of  Matrices  with  Displacement  Structure:  A  unified  treatment  of 
matrices  with  strongly  regular  and  arbitrary  rank  profiles. 

In  this  section  we  first  review  in  some  detail  our  generalized  scheme  for 
strongly-nonsingular  structured  matrices,  and  apply  it  for  illustration,  to  Hankel  and 
Quasi-Hankel  matrices.  Then  in  Section  2.2,  we  examine  the  problems  caused  by  the 
presence  of  zero  minors. 


2.1:  Factorization  of  Strongly  Nonsingular  Matrices  with  Generalized 

Displacement  Structure. 

In  1986,  Lev-Ari  and  Kailath  [1986]  extended  this  earlier  work  on  Toeplitz 
oriented  structures  to  a  much  broader  family  of  structured  matrices,  including  Hankel 
matrices  and  their  inverses,  Bezout  matrices,  sums  of  Hankel  and  Toeplitz  matrices 
and  several  other  interesting  matrices.  Motivated  by  Bistritz’s  [1984]  stability  test  for 
discrete  time  system  (which  is  a  3-term  immittance  recursion  for  inertia  computation  of 
the  related  Bezout  matrix)  Bistritz,  Lev-Ari  and  Kailath  [1987]  developed  the  so-called 
3  term  immittance  formulation  of  Levinson  and  Schur  algorithms  for  Hermitian 
Toeplitz  and  Quasi  Toeplitz  (congruent  to  Toeplitz)  matrices.  Recently  Bistritz  and 
Kailath  [1988]  developed  extensions  of  Schur  and  Levinson  algorithms  to  non- 
symmetric  Toeplitz  and  non-symmetric  Quasi-Toeplitz  matrices.  The  concept  of 
generating  functions  was  the  starting  point  for  much  of  this  analysis  and  we  shall  start 
at  that  point  here.  We  should  mention  that  all  the  results  need  the  assumption  that  all 
the  leading  minors  of  the  associated  matrices  have  to  be  nonzero  -  this  is  the 
assumption  of  “strong  nonsingularity”. 

The  key  idea  behind  all  these  fast  algorithms  is  the  concept  of  displacement  rank 
(see  Kailath,  Kung  and  Morf  [1979]),  as  extended  by  Lev-Ari  and  Kailath  [1986].  The 
concept  of  displacement  structure  and  its  properties  can  be  appropriately  described  by 
their  generating  functions.  The  generating  function  of  a  matrix  R  is  the  following 
bivariate  form,  viz. 


R( z,w)  =  [1  z  z2  ....  ]  R  [1  w  w2  ....  ]* 


(2.1) 


The  generating  function  of  a  Hermitian  matrix  with  generalized  displacement  structure 
(Lev-Ari  and  Kailath  [1986])  assumes  the  following  form, 


R  (z  ,w) 


G(z)J  G* (w) 
d(z,w) 


(2.2) 


where  J  is  any  constant  non  singular  Hermitian  matrix,  d(z,w)  is  a  Hermitian 
generating  function  like  R(z,w)  ,  viz. 
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d(z,w)  =  [1  z  z2  ....]  Jd  [1  w  w2  ....]* 


(2.3) 


where  Jd  is  a  constant  (possibly  singular)  Hermitian  matrix.  As  will  be  seen  below, 
common  choices  are  d(z,w)=l-zw*  (for  Toeplitz-related  matrices)  and 
d(z,w)  =  z-w*  (for  Hankel-related  matrices).  It  has  been  shown  in  Lev-Ari  and 
Kailath  [1986]  that  fast  factorization  of  R  is  possible  if  there  exist  matrix  functions 
©i(z)  and  constant  matrices  /,  (/  0  =  J)  such  that 


0,  (z )  A+1©i*(vv ) 


Ji  ~ 


d{z,w) 

d  (z  ,£ti)d(%i,w) 


(2.4) 


where 


Mi  =  Gfa)Ri-l&&)Gi&)  ,  (2.5) 

G(  (z )  Ji  G*(w ) 

Rt(z,w)  =  - — - i - ;  G0(z )  =  G(z)  ,  (2.6) 

d(z,w) 

and  is  arbitrary. 

They  showed  that  such  {©,  }  exist  if  and  only  if  the  {/, }  are  congruent  to  J 
and  d(z  ,w)  has  the  form 


d(z,w)  =  Y(z yf  (yv )  -  5(z )5* (w )  . 

Fast  factorization  of  R  is  obtained  via  the  recursion 

(2.7) 

where 

(*-$i)  <W*>  =  Gi(z)  &i (z ) 

(2.8) 

f  d(z  ,x- )  ] 

®i(z)  =  U  -  J  TTJl - zJiMi}  Ui  . 

[  d(z£i)d{$i,Xi)  ‘J 

(2.9) 

for  Ui 

matrices  that  satisfy 

UiJi+lU*  =  Ji 

(2.10) 

and  x,  any  complex  constant  such  that  J(t,-,t()  =  0  . 

The  choice  ^  =  0  in  (2.8)  results  in  triangular  factorization.  In  fact  (for 
%i  =  0)the  quantity  /?,■  (z  ,w )  in  equation  (2.6)  is  the  generating  function  of  the  Schur 
complement  of  the  top  left  entry  of  the  matrix  Ri_v  Other  choices  of  result  in  non 
triangular  factorizations.  It  is  crucial  for  the  existence  of  any  fast  triangular 
factorization  procedure  that  the  Schur  complement  Ri+i  inherits  the  displacement 
structure  of  /?,  .  Equations  (2.4)  -  (2.6)  represents  a  slightly  generalized  version  of  this 
fact. 


If  the  triangular  factorization  of  R  is 


R  =  LDL* 


(2.11) 


where  L  =  [I0,Ii,---,I„_i]  ,  lower  triangular  with  unit  diagonal  and  D  =  diag  (d, }  , 
then  we  can  recursively  compute  the  factorization  (2.11)  by  the  formula  (2.8) 
remembering  that 


for 


di  =  #,(0,0)  and  1,  (z )  =  z'/?,(z,0)/dj 
l/(z)  =  [  1  z  z2  ....]  1(-  . 


(2.12) 


(2.13) 


Example:  Hankei  Matrices 

Let  us  demonstrate  the  method  by  considering  the  family  of  Hankei  matrices.  As 
mentioned  earlier,  Hankei  matrices  do  possess  a  generalized  displacement  structure.  A 
Hankei  matrix  of  order  n  (  =  1,  2,  ,  ....  )  is  the  name  of  a  matrix  of  the  form 
hi+j  J/ljJo  ,  where  {  hk  }  s  are  arbitrary  (  £=0,1,  •  •  •  ,2n-2  ).  One  write  more 
explicitly: 


"„-l  = 


Since  triangular  factorization  is  nested,  factorization  of  the  leading  principal 
submatrices  does  not  depend  on  the  actual  size  of  the  matrix.  Hence  we  can  consider 
without  loss  of  generality  that  the  matrix  under  consideration  is  in  fact  semi-infinite, 
i.e.  0<i,j<oo.  We  define  its  semi-infinite  extension  Hn_l  eo  such  that  the  associated 
generating  function  H(z,w)  is  given  by, 

_  zh  (z )— w *  h*  (w) 


o 

1 _ 

h  l 

h2 

•  •  hn_ 2 

K- 

1 

hi 

h2 

h3 

•  •  K- 1 

K 

n 

l 

c 

hn-l 

K 

•  •  ^2n— 4 

h  2n- 

-3 

IA-i 

K 

hn+ 1 

•  •  h2n_2 

h2n 

-2  . 

H  (z  ,w) 


z-w 


where 


2ti-2 

h(z)=  £ hkzk 

k=Q 

Hence  H(z,w)  can  be  represented  in  the  form  (2.2)  for 

G(z)  =  [  1  zh(z)  ]  ,  J  =  antidiag  [ -j ,j ]  and  d(z,w)  =  j(z-w*)  . 


(2.14) 


(2.15) 


(2.16) 
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It  is  clear  that  d(z,w)  can  be  expressed  as  y(z )y* (w )  -  5(z )8* (w )  as  in  (2.7),  where 

Xz)  =  -^p-  and  8(z)  =  ^p-  (2.17) 

Thus  by  choosing  £  =  0  (necessary  for  triangular  factorization)  and  x  =  «>  (since 
d  (z  ,z  )=0  represents  the  real  line  here), 

G1(z)  =  G(z)0o(z)  ,  (2.18) 


where 


©oOO  = 


1  0 
z‘  *o 


,-l&n  1 


U 


0  » 


(2.19) 


where  k0  =  h0. 

Therefore  the  generating  function  of  the  Schur  complement  Hx(z,w)  of  h0  is, 

h\(z  )h  *  (yv  )-h  *  (w  )h  (z) 


H  i(z  ,w)  = 


(2.20) 


z— w 


2" -3  hk +1 

where  /ij(z)  =  jr  —  z* 


Jfc=0 


The  diagonal  element  d0  =  h0  and  the  generating 

function  l0(z)  of  the  first  column  of  the  triangular  factor  of  Hn_x  M  is  ~Z  ^ .  Therefore 

«o 


the  first  column  I0  of  the  triangular  factor  L  of  Hn_x  is  [1,  —  ,  . 

h0 


hn- 1 


]  .  It  is 


clear  from  (2.20)  that  the  Schur  complement  of  the  top  left  entry  of  a  Hankel  matrix  is 
not  Hankel  although  the  displacement  structure  of  both  is  the  same.  Matrices  with  this 
structure  are  the  so  called  Quasi-Hankel  matrices  (see  Lev-Ari  and  Kailath  [1986]).  A 
matrix  Q  with  the  generating  function 


Q(z,w)  =  /(z>g’(w)-/*(w)*(2) 


z-w 


(2.21) 


where  /(.)  and  g(.)  are  polynomials  in  z  is  defined  to  be  Quasi-Hankel.  It  can  be 
easily  shown  that  the  Schur  complement  of  the  top  left  entry  of  a  Quasi-Hankel  matrix 
is  also  Quasi-Hankel.  Thus  developing  recursions  for  factorization  of  Quasi-Hankel 
matrices  is  really  the  problem  at  hand.  The  triangular  factorization  of  a  Quasi-Hankel 
matrix  Q  with  the  generating  function  (2.21)  can  be  carried  out  by  the  recursion  (see 
Lev-Ari  and  Kailath  [1986]) 


zGul(z)  =  Gi(z)@i(z)  ;  G0(z)  =  [  f  (z)  ,  g(z)  ] 


(2.22) 


where 
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@o(z )  = 


a,1  Pi 

0  a, 


(2.23) 


where  =  sgn  [ Q ,  (0,0)]  and  , 

^  x  fi  (z  )Si  (w  )-//V  )&  (z ) 

{2i(z,w)  =  - ; -  . 

z-w 

The  coefficients  a,-  and  P(-  are  obtained  from 

[<x,  p,  ]  =  [kiQi(0,0)rmGi(0)  . 

The  formulation  (2.4)  -  (2.8)  allows  great  flexibility  in  the  selection  of  the  scalars 
£,i,ti  and  the  matrices  Ux  .  For  the  LDL*  factorization  problems  we  can  choose  {t, } 
and  {(/,  }  to  possibly  reduce  computation  (£,•  =  0  in  this  case).  As  a  matter  of  fact, 
Bistritz,  Lev-Ari  and  Kailath  [1987]  have  reduced  computational  requirements  of  the 
classical  Schur  and  Levinson  algorithms  by  choosing  proper  Ut . 


(2.24) 

(2.25) 


2.2:  Factorization  of  Structured  Matrices  with  Arbitrary  Rank  Profiles 

All  our  previous  work  on  efficient  factorization  of  matrices  with  displacement 
structure  always  assumed  certain  regularity  of  the  rank  profile  of  the  matrix.  This 
amounts  to  assuming  that  either  all  the  leading  principal  submatrices  are  full  rank  or  if 
a  leading  principal  submatrix  is  rank  deficient  then  all  the  following  leading  principal 
submatrices  are  rank  deficient  too.  However,  this  does  not  encompass  all  the 
possibilities.  Full  rank  matrices  with  singular  leading  principal  submatrices  are  often 
encountered  while  dealing  with  Hankel  and  Bezout  matrices  These  matrices  naturally 
arise  in  the  partial  realization  problem  and  the  problem  of  root  distribution  of  a 
polynomial  (see  Section  3).  Also  the  Toeplitz  matrices  associated  with  the  eigen  filter 
design  problem  and  certain  buffer  design  problems  need  not  have  regular  rank  profile. 
So  the  study  of  factorization  and  inversion  of  structured  matrices  with  arbitrary  rank 
profile  is  indeed  an  important  issue  and  as  mentioned  in  the  introduction  it  has  recently 
been  attacked  in  different  ways  by  several  researchers  (in  several  countries).  Our 
approach  is  outlined  below. 

The  basic  recursion  step  (2.8)  assumes  that  /?,•(£,•&)  *  0  .  Thus  for  triangular 
factorization  this  becomes  /?,(0,0)  =  0 <--  >  leading  principal  minors  do  not  vanish. 
The  classical  remedy  via  permutation  of  the  entries  of  would  destroy  the  underlying 
displacement  structure  of  Ri  and  hence  cannot  be  used  to  obtain  fast  algorithms. 
Vaidyanathan  and  Mitra  [1987]  have  suggested  moving  the  point  such  that 
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RifciAi)  =  0  and  the  recursion  (38)  can  be  carried  on.  However  it  would  not  be 
suitable  for  predetermined  values  of  ,  as  for  example  the  commonly  used  choice 
=  0  .  Actually  the  choice  of  4/  affects  the  computational  complexity  of  evaluating 

)  ;  the  best  choice  from  this  point  of  view  is  ^  =0  or  ^  =  ±1. 

So  we  have  chosen  a  different  path  based  on  the  fact  that  a  Hermitian  matrix  with 
singular  leading  principal  submatrices  has  a  "modified  triangular  factorization"  of  the 
form  , 

R  =  LDL*  (2.26) 

where,  L  is  lower  triangular  with  unity  diagonals  and  D  contains  non-zero  scalar  as 
well  as  block  entries  only  along  the  main  diagonal.  The  sizes  of  these  block  diagonal 
entries  can  be  determined  from  the  rank  profile  of  R  . 

So  even  if  Rt  (0,0)  is  zero,  it  may  be  possible  to  compute  GJ+](z)  from  G,  (z)  , 
where  t  is  the  size  of  the  block  in  D  whose  top  left  element  occupies  the  (i  ,i )  entry  in 

D  .  If  it  is  possible  to  do  so  then  we  shall  be  able  to  compute  the  block  Schur 

complement  of  the  txt  leading  principal  submatrix  of  /?,  if  all  the  jxj  (1  Zj  *t-l) 
leading  principal  minors  of  Rt  vanish,  directly  using  the  information  in  G,+1(z )  .  In 
fact  if  we  can  formulate  the  transition  from  G,  (z)  to  Gi+1(z)  as  a  sequence  of  t  single 
steps,  viz 

G,-(z)  — »  Gi+i(z)  — >  ....—»  Gi+1(z)  (2.27) 

then  the  process  of  computing  this  modified  triangular  factorization  can  be  kept 
completely  recursive  and  would  allow  us  to  determine  the  triangular  factor  L  ,  column 
by  column  rather  than  by  block  columns.  Truly  then  there  is  a  need  for  deriving  a 
modified  version  of  (2.8)  so  as  to  enable  us  to  compute  the  chain  (2.27) 

As  will  be  described  below,  we  have  solved  this  problem  for  Hankel  and  Quasi- 
Hankel  matrices  (see  below)  and  the  Bezout  matrices  connected  with  the  root 
distribution  problem  w.r.t  the  imaginary  axis.  We  are  currently  working  towards 
similar  results  for  Toeplitz  and  Bezout  matrices  associated  with  root  distribution 
problems  w.r.t.  the  unit  circle. 

Hankel  and  Quasi-Hankel  Matrices  with  Arbitrary  Rank  Profile 

It  is  known  that  Hankel  matrices  can  be  inverted  by  fast  algorithms  that  require 
0(n2)  operations  [see  Trench  (1965),  Berlekamp  (1968),  Massey  (1969),  Gohberg  and 
Semencul  (1972),  Kung  (1977),  Heinig  and  Rost  (1984),  Gohberg  Kailath  and 
Koltracht  (1985)  and  Citron  (1986)).  However,  the  problem  of  efficient  factorization 


of  Hankel  matrices  is  not  as  well  studied  as  the  inversion  problem.  Rissanen  (1971) 
and  more  recently  Kung  (1977)  have  considered  this  problem  in  the  context  of  partial 
realization  problems.  Very  recently  Lev-Ari  and  Kailath  (1986)  have  considered 
factorization  of  “strongly  regular”  Hankel  end  Quasi-Hankel  matrices  as  a  special 
case  of  their  work. 

Except  for  Lev-Ari  and  Kailath  1986)  other  authors  did  not  consider  the 
factorization  problems  associated  with  Quasi-Hankel  matrices.  As  far  as  factoring 
Hankel,  Quasi-Hankel  matrices  from  the  same  framework  is  considered,  Lev-Ari  and 
Kailath  (1986)  provide  the  results  directly  or  indirectly  for  matrices  that  are  strongly 
regular.  However  they  did  not  consider  matrices  that  are  not  strongly  regular. 

Rissanen  and  Kung  both  have  considered  Hankel  matrices  only.  Nevertheless  they 
have  considered  matrices  that  need  not  be  strongly  regular.  However  the  factorization 
obtained  by  these  authors  differ  from  the  form  of  factorization  (i.e.,  LDLT  where  L 
is  lower  triangular  with  unit  diagonal  and  D  is  either  diagonal  or  block  diagonal) 
considered  in  this  paper.  Rissanen  (1971)  obtains  a  LU  factorization  of  the  given 
Hankel  Matrix,  where  L  is  a  lower  triangular  matrix  with  unit  diagonal  elements  and  U 
is  an  “upper  triangular-staircase”  matrix  (see  Rissanen  (1971)  for  details).  U 
becomes  a  perfect  upper  triangular  matrix  when  the  given  Hankel  matrix  is  strongly 
regular,  russanen’s  method  is  recursive  but  it  is  not  a  fast  procedure.  Kung  (1977)  on 
the  other  hand  obtains  a  different  (HU  -  L )  factorization  of  the  Hankel  matrix  H 
where  U  is  an  upper-triangular  matrix  with  unit  diagonal  elements  and  L  is  a  lower 
“triangular-staircase”  matrix  (see  Kung  (1977)  for  more  details).  L  becomes  a 
perfect  lower  triangular  matrix  when  the  given  Hankel  matrix  is  “strongly  regular”. 
Kung’s  method  is  fast  as  well  as  recursive.  However  Kung’s  procedure  is  not 
completely  recursive  in  the  sense  that  it  makes  use  of  block  pivots  in  presence  of 
vanishing  principal  minors.  Kung’s  procedure  requires  computation  of  an  inner 
product  at  each  step.  Thus  this  procedure  is  not  completely  parallelizable. 
Determination  of  the  number  of  consecutive  zero  principal  minors  for  a  block  step  is 
not  straightforward  in  this  procedure  and  it  requires  computation  of  certain  “predicted 
Markov  parameters”  until  a  mismatch  between  these  “predicted”  and  the  “given” 
Markov  parameters  is  observed  (see  Kung  (1977)  for  details). 

Recently  Citron  (1986)  considered  the  problem  of  solving  a  Hankel  system  of 
equations.  Citron  refined  Kung’s  method  to  avoid  block  pivots  and  computation  of 
inner  products.  Determination  of  the  number  of  consecutive  zero  minors  for  block  step 
too  has  been  greatly  simplified  by  Citron  in  Citron  (1986).  However  the  problem  of 
fast  and  recursive  computation  of  a  triangular  factorization  of  the  form  LDL T  where 
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L  is  lower  triangular  with  unit  diagonal  and  D  is  either  diagonal  or  block  diagonal 
was  not  explicitly  solved  for  a  Hankel  matrix  (possibly  with  arbitrary  rank  profile). 

Below  we  describe  a  recursive  procedure  which  is  based  on  the  fact  that  the 
Schur  complement  (or  the  block  Schur  complement)  of  the  top  left  entry  (or  the  block) 
of  a  Hankel  or  a  Quasi-Hankel  matrix  is  Quasi-Hankel.  Except  for  Lev-Ari  and 
Kailath  (1986)  the  other  earlier  papers  did  not  make  use  of  this  fact.  Our  result 
generalizes  their  results  by  removing  the  requirement  of  “strong  regularity’’. 

Our  procedure  does  not  require  computation  of  inner  products  and  it  is 
completely  parallelizable.  Determination  of  the  size  of  a  block  factor  is  trivially  done 
by  counting  the  number  of  repeated  zeros  at  the  origin  of  a  polynomial  (just  as  in 
Citron(1986)). 

Since  our  approach  explicitly  computes  D  without  increasing  the  operation  count 
it  is  possible  to  utilize  our  method  to  compute  inertia  of  these  matrices  with  no  extra 

effort.  Over  and  above  since  the  block  diagonal  entries  of  D  are  either  lower  triangular 

Hankel  matrix  it  is  possible  to  compute  their  inertia  by  Iohvidov’s  (1982)  rules. 

The  recursions  (2.21)-(2.25)  for  triangular  factorization  of  the  Quasi-Hankel 
matrices  can‘t  be  used  if  QL  (0,0)  is  zero.  There  are  three  different  cases  that 
characterize  a  singularity  (i.e.  a  situation  with  Q,  (0,0)=0  ),  viz. 

(0  Qi(z,w)  =  0  .  (2.28) 

(»)  Qi(z,0)  =  0  ,  (3,(0, w)  =  0  but  Qi(z,w)*0.  (2.29) 

(Hi )  Qi  (z  ,0)  *  0  and  (2;  (0,w )  *  0  . 

The  first  case  is  trivial  and  indicates  the  end  of  the  factorization  procedure.  In  the 
second  case  the  first  row  and  the  first  column  both  are  zeros  whereas  the  matrix  Q 
isn’t  identically  zero.  This  situation  does  not  pose  a  real  problem  and  the  recursive 
procedure  (2.22)  can  be  continued  by  replacing  (2.22)-(2.23)  with  (  see  Pal  and  Kailath 
[1988]  for  more  details  ) 

*[  /,+i(0  tt+i(*>  1  =  [  fi(z)  8i(z)  ]  -  [  /,( 0)  gi(0)]  .  (2.31) 

The  corresponding  d,  =  0  and  l,(z)  can  be  arbitrarily  assigned  to  be  z‘.  The  third  case 
is  by  far  the  hardest.  We  proceed  to  analyze  this  situation  starting  with  the  following 
general  observation  about  the  generating  functions  of  Quasi-Hankel  matrices. 

“Given  a  Quasi-Hankel  matrix  Q  it  is  always  possible  to  find  two  polynomials 
/  (z )  and  g  (z )  such  that  at  least  one  of  these  two  polynomials  vanishes  at  the  origin 
while  satisfying  (2.21)’’.  The  proof  is  straight  forward.  Therefore  we  can  assume 
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without  loss  of  generality  that  /,  ( 0)  is  zero.  Then  it  can  be  easily  shown  that 
lim  =  0  too.  In  fact,  let  t  be  the  smallest  positive  integer  such  that 

Z  — »  0  Z 

lim  z~*f(z)  *  0  . 

Z  -4  o 


Then  it  turns  out  that  the  first  nonsingular  principal  submatrix  of  Qt  is  of  dimension  t . 
The  key  result  then  is  that  the  generating  function  Qcl<t (z ,w)  of  the  block  Schur 
complement  of  this  txt  block  is 


QCi,t(z  >*0 


fi+t  ( z  )8* »  (w )  ~  8i+t  (z  )ftn  (w ) 
* 

z-w 


(2.32) 


where  /t+,(z)  and  gJ+/(z)  can  be  computed  by  the  recursion 


fi+t(z)  =  z  V/ GO  , 


zgj+ii*)  =  8j(z)  ~  kjfi+t  (z)  ;  kj 


gj(  0) 

/ i  +t  (0) 


(2.33) 

(2.34) 


for  i  +t  >  j  >  i. 

Now  letting  Qc iit(z,w)  =  Qi+t (z,w)  one  gets  the  block  factorization  step 
Qi(z,w)  =  ai+,(z)a*+t(w)[  X^+t-i-y  — - — (zw*)''W  1  +  (zw*)'£)1+,(z ,w)  . 


>=o 


z  —  w 


The  generating  function  of  the  t  columns  of  the  corresponding  triangular  factor  is 
1  —  (zw  *  Y 

fi+,(z) - - - -  while  the  corresponding  block  diagonal  element  has  the  generating 

1  -  (zw  ) 

function  [  - ; - (zw  )f  1  1  ]. 

j= o  z  -  w 

A  major  application  area  for  these  results  is  to  the  problem  of  root  distribution  of 
polynomials,  which  we  address  next. 
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3:  Linear  system  stability  and  root  distribution  of  polynomials:  A  unified 
approach  via  inertia  computation  of  Bezoutians 

The  problem  of  root  distribution  of  polynomials  with  respect  to  (w.r.t.)  a  given 
curve  in  the  complex  plane  is  more  than  a  hundred  years  old.  The  two  most  important 
curves  in  engineering  application  are  the  imaginary  axis  and  the  unit  circle  as  they 
relate  to  the  stability  of  continuous  and  discrete  time  linear  time  invariant  systems 
respectively. 

Our  interest  is  in  fast  algorithms  for  solving  these  problems.  There  are  many 
such  solutions  in  the  literature,  especially  the  classical  Routh-Hurwitz  test  (Routh 
[1877])  for  the  imaginary  axis  case  and  the  Schur-Cohn  test  (Cohn  [1924])  for  the  unit 
circle  case.  Later  Marden  [1949]  and  Jury  [1964]  independently  developed  a  tabular 
form  of  the  Schur-Cohn  test.  Recently  some  new  results  have  been  reported  by 
Bistritz  [1984]  and  by  Lepschy,  Mian  and  Viaro  [1988].  Bistritz’  test  requires  only 
half  as  many  multiplications  as  for  the  Schur-Cohn  or  Marden-Jury  tests.  The  test  due 
to  Lepschy  et  al.  has  no  computational  advantage  over  the  well  known  Routh’s  test;  it 
is  just  an  alternative  test.  These  various  tests  have  all  been  derived  using  rather 
different  methods  and  perspectives.  These  derivations  do  not  provide  any  clues  as  to  a 
unified  framework,  nor  do  they  lead  us  to  the  derivation  of  equivalent  or  better  new 
tests. 

Moreover  there  are  polynomials  for  which  the  Routh  algorithm  or  the  Schur-Cohn 
algorithm  fail  to  compute  the  root  distribution,  leading  to  the  so-called  singularity 
problem;  In  fact  all  known  algorithms  for  root  distribution  suffer  from  this  problem 
and  a  thorough  discussion  of  the  singularity  has  been  lacking  for  a  long  time. 

Consider  Routh’s  algorithm  as  an  example.  Since  at  every  step  one  must  perform 
at  least  one  division  by  the  leading  element  of  the  previous  row  (previous  step),  the 
procedure  breaks  down  if  this  divisor  becomes  zero.  This  can  happen  in  two  ways: 
(i)  if  the  given  polynomial  f(s )  and  the  associated  polynomial  f(-s)  (for  real 
polynomials)  share  roots  and  (ii)  if  /  (s )  has  roots  on  both  sides  of  the  imaginary  axis 
(not  necessarily  located  symmetrically  in  pairs  w.r.t  the  imaginary  axis).  However  the 
converse  of  the  second  condition  is  not  necessarily  true.  If  all  the  elements  of  a  row 
of  Routh’s  array  become  zero,  then  condition  (i)  holds,  while  condition  (ii)  will  hold  if 
the  leading  element  of  a  row  is  zero  but  the  row  itself  does  not  vanish. 

Routh  was  aware  of  these  problems  and  he  himself  suggested  modifications  of  his 
algorithm  under  these  circumstances.  Unfortunately  Routh’s  modifications  do  not 
always  work  if  we  encounter  a  nonvanishing  row  with  a  zero  leading  element. 
Gantmakher  [1959]  suggested  a  different  way  of  handling  this  problem,  though  it  is 
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complicated  and  difficult  to  implement.  Possibly  for  the  first  time  Yeung  [1983]  came 
up  with  an  efficient  solution  for  polynomials  with  real  coefficients;  Yeung  derived  his 
solution  from  the  point  of  view  of  completing  generalized  Sturm  chains.  Later 
Delsarte,  Genin  and  Kamp  [1985]  independently  derived  the  same  solution  as  that  of 
Yeung’s  using  an  index  theory  framework.  Then  Agashe  [1985]  came  up  with  a 
complete  solution  for  both  real  and  the  polynomials,  using  Rouche’s  theorem  and  the 
principle  of  argument. 

A  very  similar  scenario  exists  in  the  unit  circle  setting.  As  a  matter  of  fact  an 
efficient  recursive  solution  for  dealing  with  a  singularity  of  the  first  kind  has  not  yet 
been  found.  Originally  Marden  [1949]  had  suggested  a  procedure  that  requires  starting 
fresh  with  a  related  polynomial  of  degree  higher  than  the  degree  of  the  polynomial 
which  poses  the  problem  of  singularity.  This  increase  in  degree  could  be  up  to  /  , 
when  the  degree  of  the  polynomial  that  poses  the  problem  of  singularity  is  either 
2/ +2  or  2t+l  .  Another  major  drawback  of  this  scheme  is  that  it  is  not  guaranteed 
that  once  a  singularity  is  bypassed,  another  one  will  not  be  encountered  in  the  process 
later.  Yeung  [1985]  suggested  another  procedure  along  the  lines  of  Marden’s  remedy. 
Yeung  reduced  the  degree  increment  from  t  to  unity  for  all  values  of  t  ;  however,  the 
second  difficulty  could  not  be  removed.  Then  in  1985  Delsarte  et.  al.  came  up  with  a 
solution  that  could  overcome  the  difficulties  in  the  work  of  Marden  [1949]  and  Yeung 
[1985].  However,  their  solution  falls  short  of  being  efficient  since  it  loses  the 
recursive  nature  in  overcoming  singularities  of  the  first  kind.  The  basis  of  Marden’s 
and  Yeung’s  derivations  was  the  principle  of  argument  where  as  Delsarte  et  al.  used 
their  recently  proposed  index  theory  of  pseudo-lossless  functions. 

It  should  be  apparent  by  now  that  there  are  very  many  different  derivations  and 
interpretations  in  both  the  regular  as  well  as  the  singular  cases.  A  recent  survey  paper 
by  Vaidyanathan  and  Mitra  [1987]  makes  a  unification  attempt.  However,  despite 
several  nice  features,  there  are  points  in  which  it  falls  short  of  being  a  definitive 
treatment.  For  example,  none  of  the  known  remedies  for  singular  cases  can  be  derived 
from  the  approach  in  their  paper.  Moreover,  Vaidyanathan  and  Mitra  [1987]  did  not 
address  the  issues  regarding  complex  polynomials. 

Recently  we  have  been  able  to  identify  a  unifying  framework  for  the  root 
distribution  procedures  (see  Lev-Ari,  Bistritz  and  Kailath  [1987])  based  on  computing 
the  inertia  of  certain  generalized  Bezout  matrices  constructed  from  the  coefficients  of 
the  given  polynomial,  generalizing  some  classical  results  of  Hermite  [1856]  (see  also 
the  fine  survey  paper  of  Krein  and  Naimark  [1937]).  The  new  point  is  that  these 
Bezout  matrices  possess  a  generalized  displacement  structure  in  the  sense  that  Kailath 


and  his  associates  have  been  studying  for  several  years.  They  have  shown  (see  e.g. 
Lev-Ari  and  Kailath  [1986]  also  see  the  previous  section)  that  displacement  structure 
can  be  exploited  to  develop  fast  procedures  (0(N2)  as  opposed  to  O  (A/3))  for 
triangular  factorization  of  such  matrices.  It  turns  out  that  application  of  these  methods 
to  the  generalized  Bezoutians  readily  yields  all  the  classical  tests  and  several 
equivalent  ones  in  a  simple  and  organized  way.  However  the  assumption  of  regularity 
(strong  regularity  of  the  rank  profile)  must  be  maintained.  At  this  point  it  will  be 
useful  to  begin  to  give  a  few  more  details. 

A  Bezout  matrix  B  admits  a  generating  function  representation  B(z,w)  given  by, 

B(z,w)  =  mOzU  -  (3.1) 

d(z,w) 


where 


B  (z  ,w  )=[1  z  z2 . ]  B  [1  w  w2 . ]* 


d(z  ,w  )=[1  z] 


a  P* 

P  5 


[1  w] 


and  the  sharp  (#)  denotes  a  suitably  defined  polynomial  transformation.  These  matrices 
B  serve  to  find  the  root  distribution  of  the  polynomial  /  (z )  w.r.t.  the  following 
partition  of  the  complex  plane, 

Q+  =  {  z;d(z,z)  >  0  } 

=  {  z;d(z,z)  =  0  }  (3.4) 

GL  =  {  z\d(z,z)  <  0  } 

More  specifically  the  inertia  of  the  Bezout  matrix  B  (i.e.  the  number  of  positive,  zero 
and  negative  eigenvalues)  indicates  how  many  roots  are  shared  by  /(z)  and  f*(z)  and 
how  many  remaining  roots  are  in  ft+  and  in  Q_  .  Q0  describes  the  curve  w.r.t.  which 
the  root  distribution  is  computed.  In  the  imaginary  axis  case,  d(z,w )  =  z+w*  and 
/#(z)  =  /  *  (-z  * ),  while  for  the  unit  circle,  d(z,w)  =  1-zw*  , 

/  (z)  =  znf*(l/z*),n  being  the  degree  of  /(z)  .  Generating  functions  of  Bezout 
matrices  are  called  Bezoutians. 

B(z,w )  in  (3.1)  can  be  compactly  represented  as, 

,  G(z)JG*(w) 

B(z,w)  =  -  ■  '  - — ~  (3.5) 

d(z,w) 

where  G  (z )  is  an  appropriate  1x2  row  vector  of  polynomials  and  J  is  an  appropriate 
2x2  non-singular  Hermitian  matrix.  As  noted  earlier  in  section-2,  this  is  in  the  so 
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called  generalized  displacement  form  of  Lev-Ari  and  Kailath  [1986],  who  showed  that 
(  see  also  the  previous  section  of  this  proposal  )  fast  triangular  factorization  of  B 
(hence  inertia  calculation  )  via  Sylvester’s  theorem  can  be  obtained  via  the  recursion 
(2.8).  The  recursion  is  started  with  G0(z)  =  G(z),J0  =  J,  and  the  updated 
Bezoutians  at  each  step  can  be  represented  as 


B{(z  ,w) 


G,(zV,G,-  *(w) 
d(z  ,w) 


In  the  above  recursion  the  parameter  is  arbitrary.  The  choice  =  0  yields  the 
familiar  triangular  (LDL*)  factorization  of  B  (as  mentioned  in  section-2),  where  as 
other  choices  yield  non  triangular  factorizations  with  the  {  5/ (£,•,£,• )  }  for  the  non 
zero  entries  of  the  diagonal  factor,  which  implies  that  the  inertia  of  B  is  easily 
computed  from  the  signs  of  the  {  Bi  (^  )  }. 

It  is  clear  that  we  have  a  great  flexibility  in  choosing  the  parameters  ,xt  and  the 
matrices  Ui  .  Specific  choices  may  help  reduce  computational  requirements  compared 
to  the  other  cases. 


For  simplicity  let  us  consider  polynomials  with  real  coefficients  for  now.  Then 
the  Bezoutian  defined  for  imaginary  axis  problems  is 

.  /(*)/•(»->  -  /(-*>/•(->  1  (36) 

Z  +  W 

Comparing  (3.5)  and  (3.6)  we  find  that 


G(z)  =  [  /  (z)  /*(-z)  ]  ,  d(z,w)  =  z+w* 


So  the  recursion  (2.8)  with  the  choice  =  0  becomes 


where 


where 


G(-+1(z)  =  G(-  (z )  0,  (z ) 


©,(z) 


f'-l 

i  +  -Ll 

l 

z  x<  J 

r  j 

(3.7a) 


(3.7b) 


Mi=G*m,-1W)Gi(0)  . 


(3.7c) 


The  choice  x,  =  j°°  (since  x,  must  be  on  the  jw-axis)  is  an  obvious  one.  Some 
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algebra  yields 


M-t  = 


Pi  Pi 
Pi  Pi 


(3.7d) 


where 


/i(0) 

2/,'(0) 


(3.7e) 


Therefore 


©;(z)  = 


1  -  p iZ 

-PiZ*1 


Piz 

1  +  p  ,-z 


-1  t/i . 


(3.7f) 


where  (/,•  can  be  any  J  unitary  matrix.  The  simple  choice  J  for  t/j  is  convenient  here 
because  with  this  choice  of  Uif  if  G,(z)  =  [  /,(z)  /,(-z)  ],  then  G,  (z  )©(z )  will  have 
the  form  [  a(z)  a(-z)  ]  ,  so  that  we  can  identify  G/+1(z)  as  [/,+i(z)  /,+ i(-z)  ]  • 


Then  the  recursions  are 


z/i+i(z)  =  /,(z)  -  PiZ_1[/,(z)  -  /i(-z)] 
z/i+i(-z)  =  PiZ_1[/i(z)  -  /i(-z)l  ~/,(-z) 


(3.8a) 

(3.8b) 


where 


Pi  = 


/i(0) 
2/  ,•  (0) 


and  the  inertia  of  B{  can  be  computed  from  the  signs  of  the  {p* }  . 

This  apparently  new  looking  recursion  is  (slightly)  different  in  form,  but 
completely  equivalent  in  the  amount  of  computation,  to  the  well  known  Routh 
recursion,  which  utilizes  the  even  and  odd  parts  of  a  given  polynomial.  The  Routh  test 
can  be  obtained  by  adding  and  subtracting  (3.8a,b).  It  can  be  directly  obtained  by 
starting  with  the  so  called  immittance-type  form  of  the  Bezoutian. 

At  this  point  it  will  be  useful  to  introduce  the  notion  of  scattering  and  immittance 
variables.  The  representation  in  (3.1)  of  B(z,w)  is  called  scattering-type,  since 


B  >0 


up  Oil's  1 

,t>.  l/(z)  I 


This  means  that  the  ratio  /#(z)//(z)  can  be  interpreted  as  a  (generalized)  scattering 
function  of  a  passive  system.  In  particular  this  ratio  is  a  continuous  time  scattering 
function  /*(-z*)//  (z)  of  a  lossless  system  when  Q0=  t*16  jw  ax*s»  an(*  a  discrete  time 
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scattering  function  rn/*(l/z*)//(z)  of  a  pseudo-lossless  system  when  Q0=  the  unit 
circle. 


Note  that  by  a  simple  (invertible)  transformation  we  can  rewrite  formula  (3.1)  as 
follows. 


B(z,w) 


a(z)b*  (w)  +  b(z)a*  (w) 
d  (z  ,w) 


(3.9) 


where 


a(z)  =  -=-[/ (z)+/#(z)]  ,  and 

b(z)  =  JL[/(z)-/‘(2)]  . 


This  representation  is  called  immittance-type,  since 


fi  50 


b(z) 


inf  -<Re- 
:ea  a(z) 


>  0 


Thus  the  ratio  b(z)/a(z)  is  positive  real  in  Q+  ,  which  for  Q.0=jw  axis  (respectively 
the  unit  circle)  corresponds  to  a  continuous-time  (respectively  discrete  time) 
impedance/admittance  (=  immittance)  function.  We  have  been  able  to  identify  the 
Routh’s  test  with  an  immittance-type  recursion  while  the  Schur-Cohn  test  with  a 
scattering  type  recursion  for  factorization  of  appropriate  Bezout  matrices. 

A  recursion  is  of  "scattering  type"  when 

Gi(z)  =  [/,(z)  ff(z)  ]  and  Jt  =  JT  =  _°J  , 

since  it  would  provide  the  scattering  type  representation  for  B,(z,w)  .  Similarly  a 
recursion  is  of  "immittance-type"  when 

G\ (z )  =  [ai(z)bi(z)  ]  and  /,•  =*//  =  [?  q]  * 

since  this  provides  the  immittance  type  representation  for  B,-  (z  ,w )  . 

In  the  case,  /#(z)=/(-z)  (considering  real  coefficients  only), 
a(z)  =  V2m(z)  and  b (z )  =  V2  n (z )  where  m (z )  and  n (z )  are  (using  standard 
notations)  respectively  the  even  and  odd  parts  of  /  (z ).  So  the  immittance  form  of 
B/(z  ,w )  is, 


Bl(ZfW)  =  2[m(z)w*(w)  4-  n(z)m*(w)] 

z+w* 


(3.10) 
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Since  we  are  only  interested  in  the  inertia  the  factor  of  2  could  be  ignored.  The 
choice  t i=joo  ,  and  £/;  =  ^  =  Jj  would  yield  for  =  0  the  following  recursion, 

z/n,+ i(z)  =  n,(z)  (3.11) 

z«l+i(z)  =  m,(z)  -  kiZ~lrii(z)  (3.12) 

where 

m(-  (0) 

=  lim  zm,-  (z  )/«,•  (z )  =  2p,-  =  - —  (3.13) 

z— >0  ^/+l(0) 

and  inertia  of  5/  can  be  computed  from  the  signs  of  {&,}  .  This  generates  the  familiar 
Routh  recursion  and  the  condition, 

{&,  }/Lq  >  0  «==>  signs  of  the  leading  <==>  stability. 

elements  of  the  Rouths’s 

array  are  all  the  same 

Consider  the  following  example.  Let  the  given  polynomial 
/(z)  =  s3  +  6s 2  +  11s  +  6  =  (s+l)(,s+2)(.s+3)  .  Then  the  algorithm  (3.11)- 
(3.12)  leads  to  the  following  computations  shown  in  array  form 


s° 

s2 

Ki 

m0(s) 

6 

6 

6/11 

m^s) 

11 

1 

121/60 

m2(s) 

60/11 

60/11 

m3(s) 

1 

The  above  is  the  Routh  array  for  the  reverse  polynomial  zn/(7z),  which  has  the 
same  root  distribution  w.r.t.  the  imaginary  axis.  Therefore  the  reversal  of  the 
positions  of  the  coefficients  is  not  material;  in  fact,  this  "reverse  Routh  array"  is  more 
convenient  than  the  original  array  of  Routh,  since  the  starting  polynomial  in  Routh’s 
original  method  depends  upon  whether  the  highest  power  of  /  (z )  is  even  or  odd,  while 
our  procedure  removes  any  such  ambiguity  at  once. 

Since  the  above  array  could  be  formed  only  using  the  coefficients  of  (.)}  ,  the 
even  polynomials,  a  three  term  recursion  is  probably  a  more  natural  form  for  Routh’s 
algorithm  rather  than  the  two  term  forms  we  have  been  describing.  The  desired  three 
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term  form  could  be  found  by  eliminating  n,(z)  from  (3. 1 1)-(3. 12), 

22ml+2<.z)  =  m,(z)  -  p ,mi+1(z)  (3.14) 

and 

m0(z)  =  m(z )  and  m^z)  =  n(z)/z  (3.15) 

Since  the  2-term  scattering  and  immittance  formulations  both  produced  the 
Routh’s  array,  so  will  the  3-term  formulations.  However  the  situation  is  not  so 
degenerate  if  is  not  on  the  jw  axis.  Keeping  real  is  necessary  to  keep  the 
successive  polynomials  real.  Then  for  a  real  £,•  *  0  ,  a  scattering  recursion  would  be 

(zHi)//+1(*>  =  //(*)  +  W-z)  (3.16) 

and 

~  (*Hi)/i+i<-*)  =  *,*/,•  (z)  +  //(-*)  (3-17) 

where 

*f  =  -fi(rb)lfi<b)  ■ 

It  turns  out  that  the  inertia  can  be  determined  from  the  signs  of  the 
{  £,-/(l  -fc,2)  }.  The  two  term  recursion  above  produces  the  following  backward 
decomposition  in  scattering  variables 

fi  ( 2 )  =  — 7T  t  (z Hi  )/«+i(z )  +  ki  (z  Hi  )/i+i(-z )  1  (3.1 8) 

1  ~ki 

This  can  be  identified  with  the  decomposition  of  Lepschy  et.  al.  [1988]  if  we  choose 
^,  =1  .  Since  only  the  indices  "i"  and  "i+l"  are  involved  in  the  computation,  Lepschy 
et  al.  in  fact  have  a  two-term  scattering  type  recursion  for  the  point  of  extraction 
^,  =  1  .  No  one  has  presented  any  test  that  would  correspond  to  a  three  term 
recursion  in  scattering  variables  for  an  arbitrary  *  0;  The  immittance-type 
recursions  are  also  not  yet  associated  with  any  known  test.  These  gaps  can  readily  be 
identified  and  filled  with  our  approach,  which  sytematically  classifies  the  vrious 
alternative  tests  (see  Table  1). 

Similar  characterization  of  tests  can  be  extended  to  the  unit  circle  problems  too. 
In  fact  the  Schur-Cohn  test  corresponds  to  a  two  term  scattering-type  recursion  and  the 
Bistritz  test  [1984]  corresponds  to  a  three  term  immittance  type  recursion.  These 
recursions  of  course  use  =  0  as  the  point  of  extraction.  A  two  term  immittance-type 
recursion  for  £,•  =  0  is  not  known  to  be  related  to  any  test  as  such.  However  such  a 
test  can  be  derived  from  the  recursions  in  Bistritz,  Lev-Ari  and  Kailath  [1987].  The 
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Imaginary-axis 


Scattering 

Immittance 

Two 

Essentially 

Essentially 

term 

Routh 

Routh 

Three 

Essentially 

term 

Routh 

Routh 
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on  the  boundary  Q0) 


Imaginary-axis 
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et.  al. 

[1988] 

Three 

? 

? 
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Unit  circle 


Scattering 

Immittance 

Two 

Schur- 

Bistritz 

term 

Cohn 

et  al.[1987] 

Three 

Szego 

Bistritz 

term 

[1939] 

[1984] 

Extraction  point  4,  =  0 
(£,■  is  not  on  the  boundary  Qq) 


Unit  circle 


Scattering 

Immittance 

Two 

? 

Lev-Ari 

term 

[1988] 

Three 

? 

Lev-Ari 

term 

[1988] 

Extraction  point  =  1 
(^,  on  the  boundary  ft0) 


Table  1  Classification  of  stability  tests 
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three  term  scattering- type  recursion  can  be  related  to  the  early  work  of  Szego  [1939] 
on  polynomials  orthogonal  on  the  unit  circle. 

Tests  involving  arbitrary  (non-zero)  extraction  point  in  the  unit  circle  context  are 
not  well  known.  Only  recently  Lev-Ari  [1988]  formulated  new  tests  which  utilize 
two-term  and  three-term  immittance-type  recursions  for  ^  =  1  .  In  a  way  this  test  is 
analogous  to  the  Routh  scheme  for  the  imaginary  axis  since  both  the  procedures 

The  above  discussion  brings  forth  the  following  categorization  of  the  tests  in 
terms  of  our  unified  factorization  framework.  This  can  be  well  described  by  the  tables 
shown  in  the  following  page,  use  an  extraction  point  on  the  separating  curve  Q0  • 

As  far  as  the  singular  cases  are  concerned,  very  little  is  known  the  about 
factorization  (inertia  extraction)  approach.  Singular  cases  correspond  to  Bezout 
matrices  with  arbitrary  rank  profile  (see  Section  2.2).  In  fact, 

Bj  (0,0)  =  0  <==>  singularity  in  Routh  recursion  . 

If  Bj.( 0,0)  =  0  ,  then  the  corresponding  leading  principal  submatrix  is  singular  and  a 
triangular  factorization  can  not  be  found  without  permutation  of  the  matrix  Bj.  . 
However  a  modified  triangular  factorization  (as  mentioned  in  section-2)  of  the  form 
Bj  =  LDL*  ,  where  L  is  lower  triangular  with  unit  diagonal  element  and  D  is  a 
matrix  with  non-zero  scalar  and  block  entries  only  along  the  main  diagonal  can  be 
computed.  The  imaginary  axis  Bezoutians  are  structurally  close  to  the  Quasi-Hankel 
matrices  (a  Bezout  matrix  is  the  product  of  a  Quasi-Hankel  matrix  and  a  signature 
matrix  (see  Pal  and  Kailath  [1988]))  it  has  been  possible  to  appropriately  modify  the 
results  in  Section  2.2  and  derive  immittance  variable  two  term  recursions  for 
overcoming  ‘singularities  of  the  first  kind’  arising  in  Routh’s  test.  This  derivation 
which  is  based  on  our  factorization  approach  produces  the  same  result  as  Yeung’s 
[1983]  test. 

The  sizes  of  these  blocks  in  D  are  determined  by  the  rank-profile  of  Bj  and 
equals  one  more  than  twice  the  number  of  "shift-rows"  in  Yeung’s  [1983]  array.  So 
the  inertia  of  Bj  is  computed  from  the  inertia  of  D  now.  The  block  entries  on  the 
diagonal  are  also  structured  and  the  inertia  of  these  blocks  can  be  just  computed  by 
formula  (see  Pal  and  Kailath  [1988]). 

The  recursion  (3.11)  -  (3.12)  gets  modified  as  follows  : 

Case  (i):  lim  m.(z)*0  ,  lim  z-1n.(z)  =  0  and 

z  — >0  z  — >0 

t  is  the  smallest  integer  such  that  lim  z“^2/+1V(z )?K) 

z  — >0 
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zmi+l  (z)  =  «,(z)  (3-19) 

zn/+i(z)  =  m;(z)  “  ^z~{2<+1)/i1(z)  (3.20) 

Case  (ii):  lim  mAz)  =  0 

z  — >0 

z/”i+i(z)  =  «,(z)  (3.2  i) 

zrt«+i(z)  =  ™«(z)  (3-22) 

These  modified  recursions  allow  us  to  determine  the  triangular  factors  and  hence 

inertia  recursively,  column  by  column  rather  than  by  block  columns. 
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4.  Fast  Parallelizable  Array  Algorithms  via  Displacement  Representations  of 
Matrices 

In  this  section  we  present  a  matrix  formulation  of  the  general  factorization 
scheme  described  in  Section  2.  By  applying  this  scheme  to  suitably  chosen  composite 
matrices,  we  obtain  in  Section  4.2  efficient  solutions  to  several  linear  algebra 
problems. 

Many  signal  progressing  problems  require  solving  large  systems  of  linear 
equations,  either  directly  or  via  (weighted)  least  squares.  The  basic  tools  for  the 
solutions  are  triangular  factorization  and  QR  factorization.  These  factorizations 
require  0(n3)  or  0(mn2)  flops  (floating  point  operations)  for  an  m  x  n  matrix,  which 
can  often  be  excessively  large.  Therefore,  attention  focuses  on  structured  matrices, 
with  an  eye  both  to  computational  reductions  and  to  implementability  in  special 
purpose  (parallel)  hardware.  Structured  matrices  arise  in  various  problems  of  linear 
time  invariant  system.  They  include  Toeplitz  and  close-to-Toeplitz  matrices,  Hankel 
and  close-to-Hankel  matrices,  and  Vandermonde  and  Krylov  matrices.  These  matrices 
have  certain  shift-invariance  properties  as  will  be  explained  shortly,  and  are  determined 
only  by  O  ( am )  parameters,  where  a  <  min  ( m  ,n ).  Therefore,  one  may  expect  that  the 
solutions  to  the  corresponding  linear  equations  should  be  obtained  with  considerably 
less  computation  than  for  general  matrices. 

Our  definition  of  structured  matrices  is  that  they  are  such  that  their  appropriately 
defined  displacement,  has  low  rank.  There  are  several  useful  definitions  of 
displacement  matrices,  but  we  shall  focus  on  two  of  them: 

V/,F^  =  A  —  F?  AFbr ,  (4.1) 

A(f/  jF^yA  —  F?  A  —  AFbT ,  (4.2) 

where  F?  and  Fb  are  chosen  matrices.  We  shall  call  the  matrices  V(F/  and 
A(F/  the  Toeplitz  and  Hankel  displacements  of  A  with  respect  to  the  displacement 
operators  {F? ,  Fb  },  respectively.  The  displacement  operators  are  chosen  so  that  the 
displacements  have  ranks  as  low  as  possible.  The  ranks  of  the  displacements 
V(f/_f‘)A  and  A^p/^yA  are  called  the  displacement  ranks  of  the  matrix  A.  The 
computational  complexity  (as  well  as  the  space  complexity)  of  fast  algorithms  is 
proportional  to  the  displacement  rank  of  the  matrix.  While  our  previous  work  focused 
mainly  on  the  Toeplitz  displacement  (4.1),  we  are  currently  studying  the  application  of 
the  Hankel  displacement  (4.2)  to  factorization  and  inversion  of  close-to-Hankel  and 
Vandermonde  matrices.  We  should  remark  her  that  Heinig  and  Rost  (1984)  have 
already  presented  some  efficient  procedures  for  the  inversion  of  such  matrices.  As 
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mentioned  in  the  introduction  their  methodology  and,  consequently  their  computational 
procedures  significantly  differ  from  ours. 

As  an  example  of  structured  matrices,  consider  a  regular  Riccati  differential 
equation, 

P(t)  =  FP(t)  +  P(t )Ft  -  P(t )HtHP  +  GGT ,  (4.3) 

where  F  e  Rnxn,  G  e  R"xl  and  H  e  Rlx" .  Then  the  stationary  solution  P  of  (4.3) 
has  a  low-rank  Hankel  displacement, 

A  (FiF)P  =  GGt  -  PHT  HP,  rank{\F  ff)  =  2. 

Similarly,  discrete-time  Riccati  equations  have  a  rank  2  Toeplitz  displacement.  We 
first  encountered  these  structured  matrices  in  the  study  of  Wiener-Hopf  integral 
equations,  which  led  us  to  the  computationally  efficient  Chandrasekhar  filters  [see, 
e.g.,  Kailath  (1972),  Kailath  (1973),  Morf,  Sidhu  and  Kailath  (1974)]. 

Another  important  example  is  the  well-known  Toeplitz  matrix.  However, 
Toeplitz  structure,  though  frequently  encountered  (stationary  process,  time-invariant 
systems,  etc)  and  well  exploited  (Levinson  and  Schur  algorithms),  is  not  invariant 
under  several  operations  arising  in  signal  processing  and  numerical  linear  algebra  - 
products,  inversion,  factorization,  Schur  complementations,  etc.  What  is  invariant  is 
the  displacement  structure.  We  were  able  to  re-derive  both  the  Levinson  and  Schur 
algorithms  using  the  low  displacement  rank  property  of  Toeplitz  matrix  T.  Notice  that 
the  displacement  V(ZjZ)T,  with  respect  to  the  {Z,Z}  has  rank  2,  where  Z  is  the 
"shift"  matrix, 


Similarly,  for  a  Hankel  matrix  H 
rank[A{ZtZ)H ]  =  2 

and  for  a  Krylov  matrix  K ,  i.e.,  the  matrix  whose  ith  row  is  c TAt~1, 
rank[A{zA')K]  =  1 

An  important  general  property  is  that  the  family  of  matrices  with  a  given 
displacement  structure  is  closed  under  inversion  and  Schur-complementation  (it  is  also 
closed  under  addition  and  in  a  slightly  extended  form,  under  multiplication). 
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Lemma  1.  For  any  nonsingular  matrix  A , 
rank[V^Fr  f>yA]  =  rank[W^Fbr 
ranJc[A(Ff  pt>y\  ]  =  rank[A^pvr  jr/ryA-1]. 


Therefore,  in  particular, 

rank[T  -  ZTZt]  =  rank[T~l  -  ZtT~1Z ]  =  2 
rank [ZH  -  HZT ]  =  rank[ZT H~l  -  H~lZ ]  =  2 


Lemma  2.  Let  the  block  matrix 


M  = 


A  B 
C  D 


have  Hankel  (Toeplitz)  displacement  rank  a  with  respect  to  lower  triangular 
displacement  operators,  {F^ ,  Fb  }.  Then  the  matrix, 


O  O 

O  S  ’ 
. 


S  —  D  —  CA~lB 


has  the  same  Hankel  ( Toeplitz )  displacement  rank  a  with  respect  to  {F? ,  Fb  }. 


Lemma  1  and  2,  whose  proof  we  omit,  suggest  that  it  may  be  possible  to 
compute  inverses  and  Schur  complements  of  structured  matrices  by  operating  on  their 
generators,  i.e.,  a  matrix  pair  X ,  Y  that  satisfies  the  displacement  equations 


^(Ff,Fb)A  —  XyT ,  or  A^f/ j?byA  —  XY^ .  (4.4) 

Doing  so  requires  O (am2)  computations,  where  a  is  the  displacement  rank  of  A, 
whereas  working  with  the  matrix  .4  itself  requires  O  (m3)  computations. 


We  recall  that  successive  Schur-complementation  can  be  used  to  produce  the 
triangular  (LU)  factorization  of  a  matrix.  We  show  how  this  fact  can  be  used  to 
obtain  efficient  algorithms  for  QR  factorization,  inversion,  regularization,  and  solution 
of  least  squares  problems.  The  key  to  obtaining  these  results  is  a  combination  of 
efficient  algorithms  for  successive  Schur-complementation  applied  to  certain 
judiciously  chosen  "composite”  (block)  matrices. 

We  shall  briefly  outline  the  resulting  so-called  Generalized  Schur  algorithms. 
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4.1  Generalized  Schur  Algorithms. 


To  efficiently  compute  the  Schur  complement  D  -  CA  lB  of  the  matrix 


M  = 


A  B 
C  D 


we  need  first  to  obtain  a  proper  generator  for  M ,  i.e.,  a  generator  of  the  form. 


O 

o 

* 

* 

o 

o 

_ i 

K 

— 

k 

*  % -  * 

X  = 

> 

S 

le 

,  Y  = 

k 

u  a 

for  the  Toeplitz  displacement,  and 


- 1 

* 

o 

O 

o 

o 

* 

- 1 

*  X 

X  -  * 

*  - - *  * 

,  Y  = 

*  X 

X  - 

X 

X 

-  X 

X 

X 

(4.5) 


(4.6) 


for  the  Hankel  displacement.  A  non-proper  generator  of  A  can  be  converted  to  a 
proper  one  in  several  ways.  One  is  applying  the  following  matrices. 


,  c2  +  s  =  1. 


(4.7) 


Proper  generators  (4.5)  or  (4.6)  can  be  obtained  by  using  the  matrix  (4.7),  or  its  special 
cases:  Givens  rotations,  hyperbolic  rotations  and  elementary  matrices.  By  post- 
multiplying  X  and  Y  with  a  sequence  of  appropriate  matrices  5,  j ,  we  can  transform  X 
and  Y  to  proper  form  with  O  (am )  computations. 

The  next  step  is  to  modify  one  column  of  X  and  Y .  Repeating  the  same  step 
(i.e.,  transformation  to  proper  form  followed  by  a  modifications  of  a  column)  r  times, 
where  r  is  the  size  of  the  square  block  A  in  M  produces  the  generator  of  the  Schur 
complement  D  -  CA~XB.  A  detailed  description  of  the  algorithm  is  provided  in  the 
Appendix. 
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4.2  Applications  of  the  Displacement  Representation  of  Composite  Matrices. 

The  above  algorithms  can  be  applied  to  the  block  matrices 

T  I  TtT  I  T\  T2  TtT  Tt  TtT  Tt 

A  =  ;  o  »  [  /  oy  t\  o  ’  L  T  1  J’  l  1  0  J’ 

to  obtain  generalized  Gohberg-Semencul  formulas  [Kailath  and  Chun  (1988)]  or, 
equivalently,  displacement  representations,  of  the  matrices, 

T~\  (TtT)-\  TlTilT2,  T(TtT)-1Tt,  (TtT)~1Tt. 

Such  composite  matrices  can  be  used  to  construct  efficient  numerical  procedures  for 
various  linear  algebra  problems.  In  the  sequel  we  describe  several  such  applications, 
involving  Toeplitz  and  close-to-Toeplitz  matrices.  Similar  applications  to  Hankel 
matrices,  upon  which  we  shall  not  elaborate,  also  yield  interesting  results  such  as  fast 
orthogonalization  of  Vandermonde  matrices. 


1.  Orthogonalization. 

Let  A  e  ROTXn  be  a  block-Toeplitz  or  a  Toeplitz-block  matrix,  and  let  us  define 
the  block  matrix, 

'  -/  A  O' 

M  =  At  O  At  .  (4.8) 

.  O  A  I 

If  we  apply  the  generalized  Schur  algorithm  to  (4.8)  after  the  mth  step,  we  shall  have 
a  generator  of 

A7 A  At 
A  I 


After  another  n  steps  of  partial  triangularization,  we  shall  have 


AT A  At 
A  I 


[R  QT  ]  + 


O  O 
O  S 


Now,  one  can  check  that  the  matrices  Q  and  R  satisfy 
A  =  QR,  QtQ=I, 


i.e.,  we  obtained  the  QR  factorization  of  A .  This  procedure  will  need  O  (a mn )  flops. 

If  one  wish  to  compute  R~]  directly,  then  one  can  perform  m  +  n  steps  of  partial 
triangularization  with  the  matrix, 
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M  = 


-I  A  O 
A7  O  I 
0  10 


because 


AT A  / 
I  O 


Rl 


[R  R't  ]  + 


O  O 
O  S 


2.  Removing  Forward  Elimination  of  Square  Systems 

If  one’s  primary  interest  in  the  factorization  is  to  solve  the  square  Toeplitz-block 
or  block-Toeplitz  system  of  equations, 


Ax  =  b. 


(4.9) 


then  one  might  want  to  obtain  the  transformed  right-side  vector  y  =  L-1b,  during  the 
course  of  the  factorization  process.  This  can  be  done  by  performing  the  following 
partial  triangular  factorization  of  the  matrix  M , 


whence  the  solution  to  (4.9)  can  be  obtained  by  solving  the  triangular  system  of 
equations, 


LT x  -  y. 


(4.10) 


3.  Removing  Back-Substitution  of  Square  Systems. 

Even  after  eliminating  the  forward  elimination  step,  from  a  hardware 
implementational  point  of  view,  the  back-substitution  step  in  (4.10)  can  still  be  quite 
cumbersome.  This  back-substitution  process  can  also  be  eliminated  by  performing  the 
partial  factorization  of  the  matrix, 


M  = 


A  -b 
I  0 


Notice  that  the  solution  x  =  T-1  b  is  the  Schur  complement  of  T  in  M .  Therefore,  after 
the  n  steps  of  partial  triangularization,  we  shall  have  a  generator  of  solution,  and  from 
the  generator,  we  can  read  out  the  solution. 

4,  Solving  Least  Squares  Problems  without  Back-substitutions. 
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To  solve  the  weighted  least  squares  problem,  that  minimizes 
p2(A,x-b)lb 

where  A  t  and  A  2  are  block-Toeplitz  or  Toeplitz-block,  we  form  the  following  matrix. 

-A  2  A  j  -b 
M  =  A\  O  0 
0  1  0 

Notice  that  the  solution, 

x  =  (A[A2_1A1rU[b  (4.11) 

is  the  Schur  complement  of 

-A  2  Ax 
A\  O 

Therefore,  after  the  m  +  n  steps  of  generalized  Schur  algorithm,  we  shall  have  a 
generator  of  the  solution  (4.11),  and  the  solution  can  be  read  out  from  the  generator. 


5.  Regularization 

If  the  given  Toeplitz  least  squares  system  is  particularly  ill-conditioned,  it  is 
meaningless  to  compute  the  exact  (least  squares)  solution,  since  small  perturbation  to 
the  matrix  can  cause  very  large  perturbations  in  the  solution.  In  such  cases,  we  solve 
the  following  regularized  system 
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APPENDIX  A 

The  following  two  theorems  describe  how  to  obtain  a  generator  of  a  desired 
Schur  complement  as  well  as  LU  factorization. 

Theorem  1.  Fast  LU  factorization  with  Toeplitz  generators. 

Let 

VjV10  =  X<»=  [*p,  x|‘>,  • .  x<»],  y<‘>=  [y[lK  yi‘>,  • .  y<»], 

where  F?  and  Fb  are  any  strictly  lower  triangular,  and  {X^,  F^}  is  proper.  Then 
the  first  column  of  L  and  the  first  row  of  U ,  of  A  ^  =  LU  is  given  by 

li  =  *i>  ui=y„ 

and  the  1st  order  Schur  complement  A  ®  -  Ijuf  has  a  generator, 

X(2)  =  [F/Xl(l),  xp,  •  •  ,  X(1)],  y(2)  =  [Fbyl'\  y?\  •  •  ,  yO)]  (A.l) 

i.e., 

V^^X*  2)y(2)7 


Theorem  2.  Fast  LU  factorization  with  Hankel  generators. 

Let 

V/^(1)  =*(1*(1)r>  *(1)=  4l)>  ■ .  4%  r(1)=  [yi(I),  yp\  • ,  y^L 

where  F?  and  Fb  are  strictly  lower  triangular  matrices,  and  {X^,  F^}  is  proper. 
Then  the  first  column  of  L  and  the  first  row  of  U,  of  =  LU  is  given  by 

Ii  =  yp(l)(F'  )Vl,(D,  ut  =  -Xi(l)(Fft)+yi, 

where  the  superscript  +  denotes  the  pseudo-inverse,  and  the  1st  order  Schur 
complement  A  2  =  A  j  -  has  a  generator, 

X<2>  =  [(x^-l,),  xf>,  •  • ,  x^)],  F<2>  =  [y P>,  yp),  •  •  ,  (y^-ti,)]  (A.2) 

i.e., 

V(f/,f‘)A(2)=X(2)F(2)7’. 


The  above  theorems  can  be  immediately  translated  into  the  following  algorithms 
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Algorithm  1. 

Input:  A  generator  of  A(1^ 

Output:  (1)  A  generator  of  the  rth  order  Schur  complement,  (2)  Triangular 

factorization. 

for  i  :=  1  to  r  do  begin 

Construct  a  proper  generator  of  At 
Obtain  a  generator  of  A(,+1)  by  (A.l); 

end 

Algorithm  2. 

Input:  A  generator  of 

Output:  (1)  A  generator  of  the  rth  order  Schur  complement,  (2)  Triangular 

factorization. 

for  i  :=  1  to  r  do  begin 

Construct  a  proper  generator  of  A 
Obtain  a  generator  of  A^'+I)  by  (A.2); 


end 
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